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An  Efficient  Reduced  Order  Modeling  Technique  for  Nonlinear 
Vibration  Analysis  of  Structures  with  Intermittent  Contact 

Akira  Saito*,  Matthew  P.  Castanier^,  and  Bogdan  1.  Epureanu^ 

Department  of  Mechanical  Engineering 
The  University  of  Michigan,  Ann  Arbor,  Ml  48109-2125,  USA 

Christophe  Pierre^ 

Faculty  of  Engineering 

McGill  University,  Montreal,  Quebec,  H3A  2K6,  Canada 

In  this  paper,  a  reduced  order  modeling  framework  for  nonlinear  vibration  problems  of  elastic 
structures  involving  intermittent  contact  is  proposed.  Of  particular  interest  is  a  vibration  problem 
of  plate-like  elastic  structures  with  a  crack  with  a  large  number  of  degrees  of  freedom  involved  on 
the  crack  surfaces.  Due  to  the  localized  nature  of  such  nonlinearity,  the  number  of  degrees  of  free¬ 
dom  on  the  surfaces  greatly  affects  the  computational  time  of  the  analysis.  Therefore,  reducing  the 
number  of  degrees  of  freedom  on  the  crack  surfaces  without  significantly  sacrificing  the  accuracy 
of  the  results  is  a  critical  issue  for  conducting  vibration  analysis  of  such  structures  in  a  reasonable 
amount  of  time.  The  focus  is  placed  on  the  development  of  an  efficient  algorithm  to  select  a  set 
of  nodes  on  the  crack  surfaces,  where  nonlinear  boundary  conditions  are  imposed.  The  method  is 
developed  based  on  a  procedure  for  selecting  master  degrees  of  freedom  for  Guyan  reduction.  The 
accuracy,  efficiency,  and  optimality  of  the  method  are  discussed  in  detail  and  compared  with  those 
aspects  of  previous  methods.  The  advantages  of  the  new  method  are  demonstrated  in  terms  of  the 
accuracy  of  the  frequency  response  and  the  resonant  frequencies. 

I.  Introduction 

VIBRATION  problems  of  structures  with  intermittent  contact  have  been  studied  extensively  for  several 
decades.  These  problems  have  practical  importance  and  feature  theoretical  complexity  due  to  their 
nonlinear  nature.  A  numerical  modeling  procedure  of  such  problems  based  on  the  finite  element  (FE) 
method  is  presented  in  this  paper.  This  work  is  motivated  by  a  need  for  developing  a  model-based  crack 
detection  algorithm  of  elastic  structures  based  on  their  spectral  properties,  such  as  resonant  frequencies  and 
response  shapes.  In  order  to  properly  predict  the  resonant  frequencies  of  such  structures,  one  has  to  consider 
the  nonlinearity  caused  by  intermittent  contact  at  the  cracks,  the  so-called  closing  crack  or  breathing  crack 
effect.  This  has  hindered  analysts  from  accurately  calculating  the  resonant  frequencies  of  cracked  structures, 
because  they  cannot  be  calculated  from  classical  linear  modal  analysis.  As  some  sophisticated  contact 
algorithms  have  been  developed,  such  as  the  penalty  method^  and  the  augmented  Lagrangian  method,^  the 
accuracy  of  the  results  of  time  transient  simulation  with  FE  models  involving  intermittent  contact  has  been 
improved.  Furthermore,  studying  vibration  problems  of  such  structures  with  an  FE  model  with  a  realistic 
complexity  is  becoming  feasible  with  the  aid  of  high-performance  computers.  However,  in  turn,  due  to 
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the  advancement  of  these  technologies,  analysts  tend  to  create  models  with  a  large  number  of  degrees  of 
freedom  (DOF).  This  is  based  on  the  expectation  that,  as  the  model  becomes  more  realistic  and  the  results 
become  more  accurate,  the  problem  can  still  be  solved  in  a  reasonable  amount  of  time.  However,  at  some 
point  the  number  of  DOF  will  overwhelm  even  the  most  advanced  hardware  and  software.  In  fact,  as  the 
model  complexity  increases,  the  cost  of  solving  contact  problems  increases  dramatically,  even  when  the 
potential  contact  areas  are  known  a  priori.  This  occurs  even  if  one  uses  reduced  order  modeling  techniques, 
such  as  the  Craig-Bampton  method.^  For  forced  response  vibration  problems  of  such  structures,  one  can 
use  accurate  and  efficient  semi-analytical  methods  such  as  the  ones  based  on  the  harmonic  balance  method 
(e.g.,  Ref.^),  by  representing  the  steady-state  dynamic  response  of  the  model  with  a  truncated  Fourier  series. 
However,  such  methods  still  suffer  from  the  increase  of  computational  cost  as  it  requires  a  fair  number  of 
harmonics  to  be  included  for  the  Fourier  transform,  in  order  to  obtain  an  accurate  result.  Therefore,  the 
goal  of  this  paper  is  to  propose  an  efficient  reduced  order  modeling  framework  for  vibration  problems  of 
elastic  structures  involving  intermittent  contact,  with  particular  attention  to  modeling  nonlinear  vibration  of 
cracked  structures.  The  focus  is  placed  upon  reducing  the  number  of  DOF  involved  in  the  contact  regions, 
in  an  automatic  manner. 

This  paper  is  organized  as  follows.  In  section  II,  a  literature  survey  over  the  related  fields  is  provided.  In 
section  III,  fhe  proposed  modeling  framework  is  presenfed,  including  fhe  reduced  order  modeling  approach 
and  confacf  DOF  selecfion  mefhod.  As  applicafions  of  fhe  mefhod,  fwo  case  sfudies  are  shown  in  section 
IV,  using  FE  models  of  a  canfilevered  cracked  plafe  and  an  academic  blade  model.  Conclusions  of  fhe  paper 
are  fhen  given  in  secfion  V. 


II.  Background 

The  issues  of  reducing  and  selecting  DOF  of  FE  models  have  been  exfensively  sfudied  by  various  mefh- 
ods  and  perspecfives,  such  as  fhe  reducfion  of  fhe  inlerface  DOE  befween  subsfruclures,  selection  of  master 
DOE  for  Guyan  reducfion,^  optimal  sensor  placemenl,  and  opfimal  consfrainf  locations.  However,  many  of 
fhe  available  mefhods  share  fo  some  exfenl  similar  goals  and  relafed  fo  each  ofher  as  described  lafer. 

Eirsfly,  fhe  issue  of  reducing  the  number  of  interface  DOF  befween  fhe  componenfs  has  been  sfudied 
by  several  researchers.  Brahmi  et  alf  proposed  a  mefhod  for  reducing  fhe  number  of  inferface  DOE  be¬ 
fore  fhe  assembly  of  subsfruclures,  where  basis  veclors  are  chosen  based  on  fhe  combination  of  secondary 
modal  analysis  of  fhe  inlerface  DOE  parlilions  of  fhe  malrices,  and  fhe  Iruncafion  of  modes  based  on  fhe 
singular  value  decomposilion.  Balmes^  inlroduced  fhe  framework  for  generalizing  inlerface  DOE  such  as 
consfrainf  modes,  by  considering  fhe  new  basis  representing  fhe  aclual  inlerface  displacemenls.  Caslanier 
et  al.  ^  also  proposed  a  lechnique  for  reducing  fhe  number  of  inlerface  DOE  by  applying  fhe  modal  analysis 
and  mode  Iruncafion  fo  fhe  consfrainf  mode  parlilion  of  fhe  malrices  produced  by  fhe  Componenl  Mode 
Synlhesis  (CMS),^  fhe  resulling  modes  of  which  are  called  fhe  characlerislic  consfrainf  modes  after  being 
Iransformed  back  info  fhe  finile  elemenl  coordinates.  All  of  Ihese  mefhods  achieve  fhe  order  reduction  of 
fhe  DOE  al  fhe  inlerface.  However,  Ihey  do  nol  provide  any  criferia  as  fo  how  fhe  inlerface  DOE  need  be 
selecfed  for  accuralely  enforcing  fhe  boundary  condifions. 

Secondly,  fhe  selection  of  master  DOF  has  been  a  crucial  factor  for  delermining  fhe  speclral  properly 
of  fhe  reduced  order  model  for  Guyan  reduclion-based  reducfion  techniques,  and  many  algorilhms  for  fhe 
selecfion  of  fhe  masler  DOE  have  been  developed.  As  if  shall  be  discussed  lafer,  Ihis  class  of  mefhods 
produces  fhe  resulls  lhal  tend  fo  solve  fhe  oplimizalion  problem  for  fhe  problems  considered  in  Ihis  sludy, 
Ihus  if  is  very  relevanl  fo  our  objeclive.  An  aulomalic  master  DOE  selecfion  algorilhm  was  firsl  proposed 
by  Henshell  and  Ong,^  in  which  fhe  masler  DOE  are  chosen  where  fhe  inertia  is  high  and  fhe  sliffness  is 
low,  whereas  fhe  slave  DOE  are  chosen  where  fhe  inertia  is  low  and  fhe  stiffness  is  high.  This  process  can 
be  aufomaled  by  examining  fhe  radian  frequency  ujg  defined  by  fixing  all  DOE  excepl  fhe  DOE  index  s. 
Namely,  ujs  =  \/kss/rnss,  for  s  =  1, . . . ,  n,  where  kij  and  rriij  are  fhe  enlries  af  fhe  ilh  row  and  jlh  column 
in  EE  stiffness  and  mass  malrices  of  size  n.  The  index  s  wilh  fhe  largesl  uig  is  idenlified  al  each  iteration 
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step,  and  the  DOF  is  eliminated  by  an  applieation  of  Guyan  reduetion  with  s  being  the  slave  DOF  and  all  the 
other  DOF  being  the  master  DOF  This  proeess  ean  be  repeated  until  the  number  of  master  DOF  reaehes  the 
desired  number.  An  approaeh  similar  to  this  algorithm  was  proposed  by  Shah  and  Raymund^°  based  on  the 
diseussion  of  Kidder  and  Flax,'^  where  the  number  of  master  DOF  is  eontrolled  by  iteratively  eliminating 
the  DOF  whose  is  larger  than  the  pre-defined  eut-off  frequeney  Wc  that  is  ehosen  to  be  approximately 
three  times  the  highest  signifieant  frequeney  in  the  frequeney  range  of  interest.  Independently  from  the  work 
by  Henshell  and  Ong,  Grinenko  and  Mokeev  developed  an  order  reduetion  teehnique  named  frequeney- 
dynamie  eondensation,^"*  whieh  also  proposed  a  eriterion  to  seleet  master  degree  of  freedom.  Although  their 
eriterion  was  legitimate,  the  implementation  of  the  seleetion  algorithm  still  suffers  from  tedious  exhaustive- 
seareh  caleulation  for  seleeting  the  DOF.  The  seleetion  method  proposed  by  Matta^^  also  uses  the  ratio 
kss/russ  with  the  similar  eriterion  to  that  proposed  by  Henshell  and  Ong.^  It  was  addressed  that  the  method 
ean  be  applied  not  only  to  the  Guyan  reduetion  but  also  to  the  CMS,  where  both  statie  and  vibration  modes 
are  used  as  basis  veetors,  onto  which  the  system  dynamics  are  projected.  A  method  proposed  by  Bouhaddi 
and  Fillod^^  used  a  different  concept  where  if  a  DOF  a  is  a  node  of  the  ith  eigenmode,  then  fixing  the  DOF  a 
results  in  Ai_i  =  Xi  where  A*  is  the  fth  eigenvalue  of  the  non-fixed  system,  and  Ai_i  is  the  eigenvalue  of  the 
system  with  the  DOF  a  being  fixed.  This  concepf  may  be  understood  using  a  vibrafion  problem  of  a  siring 
wilh  bolh  ends  fixed.  Thai  is,  Ihe  lowesl  nalural  frequency  of  Ihe  siring  wilh  a  single  supporl  becomes  Ihe 
highesl,  if  Ihe  supporl  is  placed  al  Ihe  node  of  Ihe  second  mode  of  vibralion.  This  is  because  Ihe  firsl  mode 
wilh  Ihe  conslrainl  Ihen  becomes  idenlical  to  Ihe  second  mode  of  Ihe  unconslrained  siring,  which  has  Ihe 
eigenvalue  as  Ihe  feasible  upper  bound  of  Ihe  firsl  eigenvalue  wilh  a  single  conslrainl.  Il  is  noleworlhy  lhal 
Bouhaddi  and  Fillod  explicilly  aimed  for  maximizing  the  minimum  eigenvalue  of  the  system  where  all  the 
master  DOF  fixed.  This  concepl  will  be  revisited  in  III.C.  The  melhods  for  Ihe  node  selection  reviewed  so  far 
are  based  on  Henshell  and  Ong  melhod  to  some  extent  Anolher  class  of  melhods  is  lhal  based  on  Ihe  concepl 
of  modal  energy.  The  melhod  proposed  by  Kim  and  Choi^^  uses  Ihe  energy  dislribulion  among  Ihe  DOF  for 
each  mode,  and  by  faking  Ihe  partial  sum  over  Ihe  rows  of  whal  Ihey  call  energy  dislribulion  malrix,  primary 
DOF  sel  can  be  chosen.  On  Ihe  olher  hand  Ihe  melhod  proposed  by  Cho  and  Kim^^  utilizes  Ihe  energy 
estimation  in  elemenl-level  by  Ihe  Rayleigh  quolienl  value  of  each  element  Kim  and  Cho  Ihen  proposed 
a  selection  melhod  consisting  of  Iwo  steps;^°  model  order  reduction  by  Improved  Reduced  System  (IRS)^^ 
using  Ihe  master  DOF  selected  via  a  melhod  based  on  energy  estimation  of  each  element'^  and  subsequenl 
sequential  elimination  melhod^  wilh  an  iterative  IRS.  Anolher  automatic  DOF  selection  melhod  named 
modal  energy  selection  melhod  proposed  by  Li^^  uses  melric  called  index  of  classification,  based  on  Ihe 
approximate  modal  energy  associated  wilh  each  DOF.  The  melhod  was  successfully  applied  to  an  FE  model 
of  a  cantilever  beam  problem.  Oh  and  Park^^  also  proposed  a  criterion  for  selecting  Ihe  master  DOF  based 
on  singular  values  of  Ihe  modal  malrix,  however,  il  suffers  from  Ihe  compulalional  cosl  due  to  exhaustive 
search  over  Ihe  possible  master  DOF  sels,  and  depends  on  engineer’s  knowledge  and  inluilion. 

Thirdly,  a  similar  bul  slighlly  differenl  issue  is  Ihe  selection  of  measurement  locations  for  vibration 
testing.  For  example,  one  may  need  to  measure  vibration  displacemenls  of  a  slruclure  to  determine  vibration 
modes,  typically  wilh  a  limited  number  of  sensors  or  Ihe  locations  where  Ihe  sensors  can  be  placed.  Thus, 
one  may  like  to  maximize  Ihe  information  one  can  oblain  as  much  as  possible,  wilh  Ihe  limited  number  of 
sensors  or  locations.  However  Ihe  question  arises  as  to  how  Ihe  sensors  need  to  be  located,  since  Ihe  optimal 
configuration  of  sensors  for  such  objective  cannol  be  easily  determined.  There  have  been  many  melhods 
developed  to  date  for  achieving  Ibis  goal  wilh  various  approaches.  In  particular,  one  of  Ihe  successful 
approaches  are  based  on  information  Iheory,  which  determine  Ihe  sensor  locations  by  optimizing  a  norm 
of  Ihe  Fisher  information  malrix.^"*  Among  Ihem,  one  of  Ihe  mosl  widely  used  techniques  is  Ihe  effective 
independence  vector  melhod,  or  Ihe  EIDV^^  melhod  developed  by  Kammer.^^  The  melhod  determines  Ihe 
placemenl  of  sensors  wilhin  Ihe  candidate  locations  while  mainlaining  as  much  independenl  information  as 
possible,  i.e.,  mainlaining  Ihe  measured  mode  shapes  as  independenl  as  possible.  Therefore,  il  is  nalural  to 
hypolhesize  lhal  Ihe  application  of  Ihe  nonlinear  boundary  conditions  to  Ihe  optimum  sensor  locations  would 
also  well  represenl  Ihe  real  boundary  conditions  where  Ihe  boundary  conditions  are  applied  to  all  locations 
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in  the  region.  This  method  is  hereby  eonsidered  in  this  study  and  the  formulation  is  diseussed  in  detail  in 
III.C. 

Lastly,  the  issue  of  finding  the  optimal  constraint  locations  to  maximize  the  fundamental  natural  fre¬ 
quency  of  a  strueture  is  eonsidered.  This  issue  has  an  important  relationship  with  the  optimal  master  DOF 
seleetion.  For  instanee,  suppose  there  is  a  strueture  that  ean  vibrate,  and  one  may  want  to  inerease  the  lowest 
natural  frequeney  as  mueh  as  possible,  by  alloeating  a  finite  number  of  supports  or  kinematieal  eonstraints 
to  the  strueture.  However,  the  problem  of  finding  fhe  opfimal  number  and  fhe  loeafions  of  sueh  supporfs  is 
nof  as  easy  as  if  appears.  Therefore,  if  may  be  neeessary  fo  apply  mafhemalieally  expensive  opfimizafion 
algorifhms  fo  obfain  sueh  supporf  loeafions,  sueh  as  done  in  fhe  work  by  Zhu  and  Zhang. On  fhe  ofher 
hand,  Akesson  and  Olhoff^^  sfudied  fhe  problem  by  applying  fhe  Couranf’s  maximum-minimum  prineiple. 
Namely  if  fhere  is  a  diserefe  dynamieal  system  of  size  n  and  fhere  are  r  (<  n)  kinematieal  eonsfrainfs  ap¬ 
plied  fo  fhe  system,  all  fhe  eigenvalues  of  fhe  sfruefure  inerease,  and  fhe  inereased  eigenvalues  are  bounded 
by  fhe  following  formula: 

f  G{l,2,...,n}  (1) 

where  A^  and  A°^_^  denofe  fhe  ifh  and  [i  r)lh  eigenvalues  of  fhe  sfruefure  wifhouf  fhe  eonsfrainfs,  and 
Aj  is  fhe  ifh  eigenvalue  of  fhe  eonsfrained  sfruefure.  Also  based  on  fhe  same  prineiple  and  fhe  findings  of 
Ref.,^^  Won  and  Park^*^  applied  minimizafion  mefhod  fo  obfain  fhe  opfimal  supporf  loeafion  fo  aehieve  fhe 
maximum  fundamenfal  nafural  frequeney  of  a  eanfilevered  plafe.  They  showed  fhaf  fhe  opfimal  supporf 
loeafions  should  be  on  fhe  nodal  lines  of  fhe  (r  -|-  l)lh  mode  of  fhe  uneonsfrained  sfruefure.  If  is  nofed 
fhaf  fhis  resulf  eonforms  fo  fhe  vibration  problem  of  a  fixed  siring  menfioned  above.  This  mefhod  was 
sueeessfully  applied  fo  fheir  speeifie  examples,  buf  fhe  mefhod  ean  be  applied  only  fo  speeial  eases  if  fhe 
pofenfially-eonsfrained  region  is  fhe  enfire  region  of  fhe  sfruefure,  where  fhe  poinfs  on  fhe  nodal  lines  ean 
be  seleefed.  Namely,  if  fhe  regions  fo  whieh  fhe  eonsfrainfs  are  applied  are  limifed  fo  some  speeifie  regions 
of  fhe  sfruefure,  fhen  fhe  nodal  lines  may  nof  exisf  in  sueh  regions,  and  fhe  minimizafion  problem  beeomes 
more  eomplieafed. 

If  is  inferesfing  fo  nofe  fhaf  fhe  idea  of  eonsfraining  fhe  nodal  lines  was  used  fo  opfimally  seleef  fhe  master 
DOF  for  Guyan  reduefion  by  Bouhaddi  and  Fillod,^^  buf  fhey  were  nof  aware  of  fhe  applieabilify  of  fheir 
mefhod  fo  opfimally  seleef  fhe  supporf  positions,  while  Won  and  Park  were  nof  aware  of  fhe  applieabilify 
of  fheir  mefhod  fo  opfimally  seleef  fhe  masfer  DOF  loeafions  for  Guyan  reduefion.  In  Ibis  paper,  we  lake 
advanlage  of  Ibis  similarily  belween  fhe  optimal  master  DOF  seleefions  and  fhe  eonsfrainf  loeafions,  in  order 
fo  aehieve  fhe  opfimal  seleetion  of  fhe  DOF  where  fhe  nonlinear  boundary  eondifions  are  applied. 

III.  Mathematical  Formulation 

Consider  small  vibrafion  problems  of  an  elasfie  sfruefure  represented  as  O  wilh  a  fixed  boundary  F^, 
where  fhe  sfruefure  may  or  may  nof  involve  intermittent  eontaet  at  F ^  and  F b  during  the  vibration  eyeles, 
sueh  as  shown  in  Fig.  1.  Namely  the  boundaries  open  and  elose,  thus  the  vibration  problem  is  nonlinear  be- 
eause  the  eondition  for  the  the  boundaries  to  be  in  eontaet  is  dependent  on  the  displaeement  field  ilself.  Thai 
is,  fhe  boundary  eondifions  al  Fyi  and  F^  are  nonlinear.  If  is  well  known  fhaf  fhe  syslem  eigenveefors  and 
eigen-frequeneies  are  differenl  from  fhe  aelual  response  shapes  and  resonanl  frequeneies  of  Ibis  nonlinear 
problem.  In  Ibis  paper,  fhey  are  respeelively  referred  fo  as  fhe  nonlinear  normal  modes  (NNMs)  and  NNM 
frequeneies,  as  was  also  done  in  Ref.^^ 

Now,  if  fhe  sfruefure  is  diserelized  wilh  a  mefhod  sueh  as  finife  elemenl  mefhod,  fhe  nonlinearily  asso- 
eialed  wilh  fhe  eonlael  is  localized,  in  fhe  sense  fhaf  fhe  nonlinearily  is  eaused  only  by  a  small  porlion  of 
fhe  enfire  sfruefure.  In  fhe  following  formulalions,  a  sel  of  indiees  of  DOF  in  sueh  region  is  denoted  as  B 
(boundary),  whereas  a  sel  of  indiees  of  fhe  DOF  in  fhe  resl  of  fhe  regions  is  denoled  as  X  (infernal),  and 
partitions  of  veelors  and  malriees  assoeialed  wilh  Ihese  sels  are  designaled  wilh  subseripls  of  fhe  assoeialed 
lower-case  ilalic  tellers,  i.e.,  b  and  i.  The  sizes  of  fhe  sels  are  denoled  as  \B\  =  ng  and  |X|  =  nj.  All  fhe 
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Figure  1.  An  elastic  structure  with  potentially  contacting  boundaries 


Other  DOF  sets  defined  hereinafter  follow  the  same  notation. 

If  the  finite  element  mass  and  stiffness  matrices  are  denoted  as  M  G  and  K  G  and  the  nodal 
displacement  vector  is  given  as  x  G  the  governing  equations  of  the  vibration  problem  with  the  absence 
of  external  forcing  and  damping  may  be  written  in  a  partitioned  matrix- vector  form  as  follows: 


f  Mftfe  MfeA 

^  / Kbb  KfeA 

Xft 

ffeCxfc) 

{Mib  MiiJ 

Xj 

\Kib  KiiJ 

Xj 

0 

where  a  dot  ( ’ )  denotes  a  time  derivative,  and  ff,  G  M”®  denotes  the  nonlinear  force  associated  with  the 
intermittent  contact.  When  dealing  with  this  type  of  nonlinear  vibration  problems,  one  can  apply  linear 
reduced  order  modeling  techniques,  such  as  Guyan  reduction,^  system  equivalent  reduction  expansion  pro¬ 
cess  (SEREP),^^  iterated  improved  reduced  system  (IIRS),^^’^^  or  Component  mode  synthesis  (CMS).^  With 
such  methods,  one  can  obtain  smaller  system  matrices  by  reducing  the  size  of  Xj  by  means  of  Rayleigh-Ritz 
coordinate  transformation  comprising  of  various  basis  vectors  such  as  static  deformations  and  vibration 
modes,  yet  keeping  the  accessibility  to  the  physical  coordinates  of  x;,.  Eor  instance,  with  the  help  of  CMS, 
one  can  obtain  a  system  with  desired  spectral  properties  and  accessibility  to  x^,  the  size  of  which  is  as 
small  as  ng  DOE  plus  the  number  of  linear  normal  modes  whose  frequencies  lie  in  the  frequency  ranges  of 
interest.  The  use  of  such  linear  reduced  order  modeling  methods  greatly  helps  ones  to  analyze  the  dynamic 
response  of  systems  with  localized  nonlinearities,  such  as  transient  dynamic  analysis,^^  and  nonlinear  har¬ 
monic  response  analysis. However,  even  with  these  reduced  order  modeling  methods,  if  the  number  of 
DOE  involved  in  the  b  partition  becomes  large,  especially  the  cases  with  very  fine  mesh  in  the  contacting 
regions,  one  cannot  take  advantage  of  the  linear  reduced  order  modeling  techniques,  as  the  computational 
cost  associated  with  the  nonlinear  dynamic  analysis  typically  grows  as  the  number  of  DOE  in  the  b  partition 
increases.  Eurthermore,  if  one  simply  attempts  to  eliminate  some  of  the  DOE  in  the  b  partition,  it  results 
in  inaccurate,  or  even  wrong  results,  in  comparison  to  the  results  obtained  with  a  full  set  of  DOE  in  the  b 
partition.  Therefore,  in  order  to  obtain  accurate  computational  results,  such  as  those  of  nonlinear  forced 
response,  one  needs  to  keep  as  many  boundary  DOE  as  possible,  which  could  easily  result  in  prohibitively 
costly  calculations.  Typically  as  a  “workaround”  to  avoid  the  inaccurate  results  due  to  the  lack  of  sufficient 
DOE  considered  and  at  the  same  time  to  obtain  efficient  computational  model,  one  has  to  select  the  DOE  in 
a  heuristic  way,  which  greatly  depends  on  the  system  characteristics  and  analyst’s  experience  and  intuition. 
Moreover,  if  the  model  is  developed  in  such  ways,  the  error  contained  in  the  following  analysis  results  can¬ 
not  be  estimated  a  priori.  Our  aim  here  is  to  develop  an  automatic  way  to  select  the  DOE  in  B  for  a  desired 
number  of  DOE  to  be  selected. 


5  of  21 


American  Institute  of  Aeronautics  and  Astronautics 


III.A.  Primary  Model  Reduction 


In  order  to  reduce  the  number  of  DOF  included  in  X  to  make  the  subsequent  development  more  efficient,  first 
a  model  reduction  is  applied  to  Eq.  (2).  Namely,  X  is  further  divided  into  two  sets,  i.e.,  T  =  O  U  D  where  O 
is  a  set  of  DOF  indices  associated  with  the  nodes  to  be  used  in  the  following  analysis,  such  as  observing  the 
behavior  of  the  system  or  applying  external  loading,  and  X)  is  the  rest  of  DOF  indices  in  X,  which  is  to  be 
apparently  deleted  from  the  system  by  the  reduction  methods.  In  addition,  a  set  of  DOF  indices  to  be  used 
as  the  master  DOF  is  defined  as  active  DOF,  designafed  as  A,  and  A  =  B  U  O.  Now  consider  an  eigenvalue 
problem  of  fhe  sysfem  Eq.  (2),  where  fhe  eigenvalue  A  and  fhe  corresponding  eigenvector  c/)  musf  satisfy  fhe 
following: 


Kaa 

0a 

=  A  ( 

0a 

Kda 

KdJ 

4>d 

\Mda 

Mm 

4>d 

(3) 


where  tp  =  In  fhis  sfudy,  a  mixed-boundary  CMS  of  Hinfz-Herling^^’^^  is  chosen  for  fhe 

primary  model  reduction.  Namely,  wifhouf  fhe  presence  of  rigid  body  modes,  fhe  coordinafe  Iransformafion 
is  defined  as 


Xa 

Xd 


Urj 


Va 

'Hm 


(4) 


where  =  rja,  Tfm  is  a  vector  of  modal  coordinafes,  'I'  and  $  are  so-called  consfrainf  modes  and  fruncafed 
free-inferface  normal  modes  in  a  modified  form,  which  are  respecfively  defined  as 


^  = 


$  = 


/ 

0 


(5) 

(6) 


and  $  =  [0(1),  0(2))  •  •  • )  0(A:)])  k  <  n,  each  subscripf  in  parenfheses  denoting  fhe  corresponding  mode 
number.  Using  fhe  fransformafion  defined  as  Eq.  (4),  fhe  projecfed  eigenvalue  problem  is  obfained  as 


KhV  =  (7) 

where  and  K.h  =  H^KH.  If  should  be  noted  fhaf  fhe  projecfed  eigenvalue  problem 

Eq.  (7)  possesses  af  teas!  fhe  same  eigenvalues  of  fhe  original  systems,  i.e.,  A(i),  A(2), . . .  X(k)  (the  indices 
may  be  differenl  from  fhe  ones  for  fhe  projecfed  sysfem.)  This  is  because  fhe  subspace  spanned  by  fhe 
columns  of  confains  fhe  eigenvectors  of  Eq.  (3),  i.e.,  0(j)  G  span('I',$)  for  j  =  l,.../c,  as 

span(^',  $)  =  span(^',  $),  and  hence  fhe  projecfed  eigenvalue  problem  has  fhe  same  eigenvalues  as  fhe 
original  ones.  If  means  fhaf,  fhe  eigenvalues  of  fhe  projecfed  sysfem  Eq.  (7)  does  not  confain  any  error  in  fhe 
eigenvalues  and  eigenvectors,  wifh  respecf  to  fhose  of  fhe  original  eigenvalue  problem  of  fhe  finife  elemenf. 
Alfhough  fhis  advanfage  comes  wifh  fhe  expense  of  calculafing  fhe  eigenvalues  and  eigenvecfors  of  fhe  finite 
elemenf  model,  if  is  nol  a  major  drawback  considering  fhaf  fhe  compufafional  cosf  involved  in  fhe  nonlinear 
compufafions  wifh  fhe  original  finife  elemenf  would  be  more  prohibifively  expensive,  fhan  calculafing  a  few 
normal  modes  of  fhe  finife  elemenf  model. 


III.B,  Nonlinear  DOF  sampling 

Wifh  fhe  reduced  order  model  obfained  in  III.A,  fhe  nexf  step  is  to  selecf  fhe  DOF  in  B  such  fhaf  fhe  nonlinear 
characferisfics  of  fhe  sysfem  can  be  well  approximafed  by  applying  fhe  nonlinear  boundary  condifions  only 
on  fhe  selecfed  DOF. 
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Figure  2.  Schematic  of  the  node  sampling:  •,  selected  node  (M) 


As  mentioned,  accurately  calculating  the  NNM  frequencies  is  of  primary  interest  of  this  study.  The  NNM 
frequencies  of  the  system  can  be  obtained  in  several  ways,  such  as  time  integration  of  Eq.  (2)  for  harmonic 
loading,  or  harmonic-balance-based  frequency/time  domain  analysis.^®’ It  was  shown  by  the  authors  that 
the  NNM  frequencies  for  cracked  plates  obtained  by  the  nonlinear  harmonic  response  analysis  can  be  well 
approximated  by  the  application  of  bilinear  frequency  approximation  even  when  the  crack  surfaces  involve 
multiple  DOF.^°  Therefore,  as  a  measure  to  evaluate  the  results  obtained  with  the  selected  DOF,  bilinear 
frequency  is  used  in  the  following  development.  Namely,  the  ith  NNM  frequency  iOni  can  be  approximated 
by  a  bilinear  frequency  oju  defined  as 


‘IuJqUJs 

LUhj  — 

1^0  + 


(8) 


where  uJo  and  ujg  are  the  natural  frequencies  of  the  corresponding  linear  systems,  which  can  be  respectively 
obtained  by  solving  the  following  eigenvalue  problems: 


Kh?7  =  (a;o)M//r7,  subject  to  open  B.C.’s  (9) 

subject  to  B.C.’s  (10) 


The  open  B.C.  is  a  boundary  condition  where  no  constraint  is  imposed  on  the  nodes  on  Tyi  and  T^,  or  the 
DOF  in  B.  Thus  in  fact  Eq.  (9)  is  identical  to  Eq.  (7).  On  the  other  hand  for  the  sliding  B.C.,  it  is  assumed 
that  Tyi  can  freely  slide  with  respect  to  T^  but  cannot  separate  along  the  normal  directions,  as  described  as 
follows. 

Here  a  contact  pair  is  defined  as  a  pair  of  nodes  on  T^i  and  T b,  which  may  or  may  nof  be  in  confact 
during  fhe  vibration,  and  a  sef  of  numbers  denofing  all  fhe  confacf  pairs  is  defined  as  Cep-  For  fhe  jfh 
confacf  pair  in  Cep,  fhree  mufually  perpendicular  normal  vectors  af  a  node  on  T^r  are  defined  as  n{,  n^, 
and  rig  where  is  fhe  normal  veefor  poinfing  oufward  from  fhe  surface,  and  rig  are  unif  vectors  fhaf 
are  fangenf  fo  fhe  surface  and  perpendicular  fo  each  ofher.  Using  fhese  veefors,  a  coordinafe  fransformafion 
mafrix  P-^  =  (rig,  n^,  rig)  is  defined  for  each  confacf  pair,  wifh  fhe  assumpfion  fhaf  a  nodal  displacemenf 
veefor  confains  only  franslafional  DOF,  such  fhaf  fhe  xi  componenf  of  fhe  displacemenf  veefor  of  fhe  node 
is  aligned  wifh  nj,  and  pointing  outward  from  the  surface.  For  the  other  node  of  the  jth  contact  pair  on  T^, 
the  corresponding  coordinate  transformation  matrix  that  aligns  the  xi  component  of  the  nodal  displacement 
vector  with  the  normal  vector  is  defined  as  P'^  =  — P'^.  Now  assembling  P-^  and  P-^  for  all  j  €  Cep,  a 
coordinafe  fransformafion  is  defined  as 


P  = 


^P/j  0  0  \  nc^p 

0  lo  0  ,  where  P;,  =  A(P^,P^)  (11) 

I  J=i 

V  0  0  1^/ 


and  A  is  an  assembly  operator,  Pfo  G  g  and  Im  £  Nexf,  for  fhe  jfh  confacf 

pair,  fhe  xi  componenfs  of  fhe  nodal  displacemenf  veefors,  which  are  denofed  as  and  rf^,  are  fransformed 
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to  a  relative  displaeement  =  {r]\  +  '>1b)/ ^  ^  displaeement  =  {rf^  —  rf^) /\/2.  Namely,  denoting 

the  set  of  DOF  eorresponding  to  '^A  and  for  all  yth  eontaet  pair, 


- 1 

1 - 

1 _ 

_ 1 

x/2(l  -l) 

(12) 


Now  defining  sets  X,  y,  and  Z{X[jyuZ  =  B  and  nx  =  ny  =  nz  =  nc^p)  that  respeetively  eontain  sets 
of  indiees  of  the  DOF  eorresponding  to  xi,  X2,  and  X3  for  all  j  G  Cep,  and  denoting  the  eoeffieient  matrix  in 
the  Eq.  (12)  as  R-^  ,  one  ean  define  a  fransformafion  mafrix  R  by  assembling  R-^  for  j  G  Cep  as  follows 


C-cp 

where  R^  =  A  (R^)  (13) 

f=i 

and  Ra;  G  W^xxnx _  Considering  fhaf  and  R~^  =  R"*",  fhe  eigenvalue  problem  Eq.  (7)  ean  be 

Iransformed  fo 

Kq  =  (t(;Q)lMq  (14) 

where  r/  =  PRq,  M  =  (PR)"'"Mj:fPR,  and  K  =  (PR)^K//PR.  Now  nofing  fhaf  q  ean  be  parfifioned 
info  q  =  [q^ ,  q^]  where  q^  is  fhe  veefor  of  relafive  DOE,  or  ,  Vy  G  Cep,  and  q^  is  defined  fo  be  fhe 
generalized  infernal  DOE  eonfaining  ,  Vj  G  Cep,  X2  and  X3  eomponenfs  of  fhe  nodal  displaeemenf  veefors 
of  fhe  nodes  of  fhe  eonfaef  pairs,  displaeemenf  veefors  of  fhe  observer  nodes,  and  modal  eoordinafes.  Thai 
is,  Eq.  (14)  ean  be  wrillen  as 


qr 

=  (“^o) 

qr 

^gg) 

5a 

\Mgr 

^gg) 

5a 

(15) 


where  TZ  C  X  and  Q  =  (A\7^)  U  Ad. 

The  besf  approximation  fo  fhe  NNM  frequeney  ean  be  obfained  when  fhe  sliding  boundary  eondifions 
are  imposed  on  all  of  fhe  nodes  on  fhe  surfaees  F^i  and  F b-  Namely  fhe  assoeiafed  eigenvalue  problem  wifh 
fhe  sliding  boundary  eondifions  ean  be  obfained  by  eonslraining  all  fhe  relafive  DOF,  or  q^  =  0,  i.e., 

^gg^g  =  i‘^s)^gg(ig  (16) 


Now,  we  assume  fhaf  we  do  not  like  fo  eonsider  all  nodes  in  TZ  for  fhe  subsequenf  foreed  response  analysis 
due  fo  fhe  large  number  of  DOF  involved  in  TZ.  In  ofher  words,  fhe  nodes  where  fhe  nonlinear  boundary 
eondifions  are  applied  should  be  sampled  sueh  as  illusfrafed  in  Fig.  2.  The  seleefed  DOF  is  designafed  as 
nonlinear  DOF,  and  a  sef  of  indiees  of  fhe  nonlinear  DOF  is  denofed  as  M,  where  M  dTZ.  The  resf  of  DOF 
in  TZ  is  designafed  as  linear  DOF,  and  assoeiafed  sef  is  denoted  as  C  where  M  VJ  C  =  TZ.  Therefore,  fhe 
bilinear  frequeney  should  be  ealeulafed  wifh  ujs  sueh  fhaf  fhe  sliding  B.C.  is  applied  only  on  fhe  DOF  in  M, 
or  q„  =  0,  i.e.. 


K3Z  Kgg 


qz 

^g 


qz 

^g 


(17) 


Considering  fhaf  fhe  value  of  fhe  nafural  frequeney  of  fhe  sysfem  wifh  fhe  open  boundary  eondifions,  cOq, 
is  independenf  on  neifher  fhe  number  nor  fhe  pattern  of  fhe  seleefed  DOF  (reealling  fhaf  span('J'c,  con- 
fains  fhe  ehosen  eigenveefors),  one  ean  see  from  Eq.  (8)  fhaf  tObi  is  dependenf  only  on  cOg  for  a  fixed  lVq. 
Now  eonsidering  fhe  Rayleigh’s  fheorem  of  eonsfrainfs  defined  by  Eq.  (1),  if  is  known  fhaf  all  fhe  sysfem’s 
eigenvalues  inerease  if  a  single  eonsfrainf  is  imposed  on  a  system.  Therefore,  as  fhe  number  of  eonsfrainfs 
on  Eq.  (7)  fo  ealeulafe  uos  inereases,  lOs  inereases.  Eurfhermore,  eonsidering  fhaf  LObi  is  a  monofonieally 
inereasing  funefion  of  cOg  for  a  fixed  tUo,  or  dcobildcog  =  2i0gl{u;o  +  w^)^  ^  0,  one  ean  sfafe  fhaf  fhe  besf 
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approximation  of  uibi  for  a  given  number  of  n_sf  can  be  obtained  when  the  maximum  cvg  is  achieved.  Thus  a 
corresponding  maximization  problem  is  stated  as  follows: 

max  ujs  (Af) 

Men  (18) 

subject  to  |AA|  =  nj\j- 

This  maximization  problem  may  be  solved  by  mathematical  programming  methods,  such  as  integer  pro¬ 
gramming  or  topology  optimization  methods  as  was  done  in  Ref.^^  As  it  shall  be  discussed  next,  this 
maximization  problem  can  in  fact  be  treated  in  a  more  efficient  way  by  the  use  of  Guyan  reduction  and 
some  methods  to  choose  the  master  DOF  for  reduced  order  modeling  techniques. 


III.C.  Automatic  master  DOF  selection 


The  methods  for  automatically  selecting  the  master  DOF  for  the  Guyan  reduction  have  been  previously 
developed,  such  as  in  Refs.^’  In  particular,  the  method  proposed  by  Henshell  and  Ong^  appears  to 
be  the  most  successful  approach.  Although  it  has  been  known  to  be  computationally  expensive  due  to  the 
nature  of  eliminating  a  single  DOF  per  iteration  and  successive  application  of  Guyan  reduction,  this  can 
be  alleviated  by  the  application  of  the  primary  model  reduction  by  the  CMS  as  developed  in  UFA.  As  was 
mentioned  by  Bouhaddi  and  Fillod,^^  and  Shah  and  Raymund,^°  the  master  DOF  of  Guyan  reduction  should 
be  chosen  such  that  the  valid  eigenvalue  range  of  the  reduced  order  model  is  maximized.  In  general,  it  has 
been  known  that  the  eigenvalue  range  of  validity  is  “bounded”  by  the  lowest  eigenvalue  of  the  system  with 
all  the  master  DOF  fixed.  Here  this  concept  is  applied  to  the  problem  of  finding  the  optimal  M  that  solves 
Eq.  (18).  Namely,  the  corresponding  eigenvalue  problem  is  Eq.  (17)  by  regarding  as  the  master  DOE. 
As  was  discussed  in  the  Ref.,^^  the  error  bounds  in  the  hh  eigenvalue  of  the  reduced  model  produced  by  the 
Guyan  reduction  can  be  obtained  a  priori  by  the  following  relationship 


0  ^ 


A* 


A 


s,min 


\i 


(19) 


where  =  (A*  —  A*)  / A*  is  the  relative  error  in  the  ith  eigenvalue.  A*  is  the  hh  eigenvalue  of  the  reduced  order 
model,  Aj  is  the  ith  eigenvalue  of  the  original  finite  element  model,  and  Xs,min  is  the  smallest  eigenvalue  of 
the  system  with  all  the  master  DOF  fixed.  Eor  Xi/Xs,min  ^  the  upper  bound  asymptotically  converges  to 
the  following  value,^^ 

0  ^  ^  Xi/Xg^^in  (^fl) 


Therefore,  it  is  apparent  that  maximizing  Xg^min  results  in  minimizing  the  upper  bound  of  the  error  for  all 
the  eigenvalues  of  the  reduced  order  model.  Hence  this  gives  us  a  guideline  for  selecting  the  master  DOE  for 
Guyan  reduction  such  that  the  errors  in  the  eigenvalues  of  the  resulting  reduced  order  model  are  minimized. 

By  observing  this  fact  from  another  point  of  view,  one  may  see  that  if  a  certain  set  of  master  DOE 
can  achieve  the  maximum  Xg^min,  we  can  obtain  not  only  an  accurate  reduced  order  model  that  can  well 
approximate  the  first  few  lowest  eigenvalues  of  the  original  system,  but  also  as  a  “byproduct”,  a  good 
estimate  on  the  optimal  constraint  locations  that  maximize  the  fundamental  frequency.  Recasting  this  to  our 
original  problem  of  selecting  the  optimal  set  M,  the  error  bounds  Eq.  (20)  associated  with  the  eigenvalue 
problem  Eq.  (17)  are  written  as 


0  ^ 


(21) 


where  e*  =  is  the  fth  eigenvalue  of  a  reduced  order  model,  (a;s)i  is  the  lowest 

natural  frequency  of  Eq.  (17).  The  corresponding  maximization  problem  is  Eq.  (18),  and  by  solving  this 
problem  for  the  lowest  eigenvalue,  (wsji,  one  can  expect  that  the  chosen  nodes  pattern  is  at  least  sub- 
optimal. 
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Algorithm  1  DOF  selection  based  on  Henshell  and  Ong  method 

1:  for  i  =  1  to  i  =  rj-T^  —  do 

2:  Calculate  ^ kjj / mjj  for  j  £TZ 

3:  Find  qi  such  that  \/ kq^q^  / ruq^q^  =  max 

j&TZ 

4:  £  <—  {qi, . . . ,  qn^}  where  g’s  are  the  DOF  associated  with  the  fcth  contact  pair  {k  G  Cep)  and  is 

the  number  of  DOF  in  /cth  contact  pair 
5:  M  ^  n\c 

6: 

7:  Calculate  constraint  modes  Eq.  (22) 

8:  Apply  Guyan  reduction  to  the  system  matrices:  M  ^  K  ^ 

9:  end  for 

According  to  Refs.,^^’^^  the  sequential  elimination  method  by  Henshell  and  Ong^  tends  to  keep  As, mm 
high,  as  it  eliminates  the  DOF  associated  with  the  highest  constrained  frequency  at  each  iteration  as  the 
slave  DOF.  Namely  after  the  elimination  procedure,  if  the  chosen  master  DOF  are  all  fixed,  the  system  is 
left  with  the  DOF  that  were  chosen  as  the  slave  DOF  that  were  identified  fo  have  fhe  highesf  consfrained 
frequency  af  each  elimination  process.  Thus  fhe  resulfing  sysfem  wifh  all  fhe  masfer  DOF  fixed  fends  fo  have 
a  larger  Xs,min  than  fhaf  calculafed  wifh  sysfems  wifh  ofher  possible  combinafions  of  masfer  DOF  fixed.  The 
Henshell  and  Ong’s  mefhod  fhaf  is  adapfed  specifically  for  fhis  problem  is  shown  in  Algorifhm  1.  Firsf,  af 
each  iferafion,  fhe  rafios  of  fhe  diagonal  ferms  of  fhe  sfiffness  mafrix  kjj  fo  fhe  diagonal  ferms  of  fhe  mass 
mafrix  rrijj  are  calculafed  for  Vj  G  TZ,  and  fhe  index  qi  fhaf  gives  fhe  maximum  rafio  among  j  G  72.  is 
obfained.  Nexf,  fhe  sef  C  is  updafed  such  fhaf  if  confains  qi,  and  all  fhe  ofher  DOF  fhaf  are  associafed  wifh 
fhe  confacf  pair  k  G  Cep  to  which  the  (71th  DOF  belongs,  e.g.,  the  DOF  that  are  perpendicular  to  the  normal 
direction.  The  set  M  is  then  updated  such  that  it  excludes  the  selected  DOF  of  C  from  72,  and  the  set  72  is 
re-defined  as  M.  A  consfrainf  mode  is  calculafed  by  solving  a  problem  where  a  unif  displacemenf  is  applied 
fo  a  DOF  in  M  whereas  all  fhe  ofher  DOF  in  M  being  fixed.  This  is  repeafed  for  all  DOF  in  M,  resulfing  in 
fhe  following  mafrix: 

-4  \ 

where  is  fhe  mafrix  of  consfrainf  modes  for  all  DOF  in  M.  The  Guyan  reduefion  is  fhen  applied  fo  fhe 
mass  and  sfiffness  mafrices.  The  iferafion  confinues  until  fhe  number  of  DOF  in  M  reaches  fhe  specified 
value  of  njq-  using  'J'. 

In  order  fo  clarify  fhe  appropriafeness  of  fhe  algorifhm  in  Algorifhm  1  fo  fhis  problem,  anofher  algorifhm 
for  selecting  DOF  is  shown  here  for  comparison.  The  mefhod  of  effeefive  independence  veefor,  or  fhe  EIDV 
mefhod  developed  by  Kammer,^^  is  a  mefhod  fo  choose  fhe  sensor  placemenf  locations  for  fhe  vibration 
measuremenf  of  large  scale  sfruefures.  The  mefhod  aims  fo  make  fhe  measured  eigenveefors  as  linearly 
independenf  as  possible.  According  fo  Penny  et  many  of  fhe  criferia  for  choosing  fhe  masfer  DOE 
for  model  order  reduefion  have  similar  criferia  for  choosing  measuremenf  locations  in  a  way  such  fhaf 
fhe  lower  frequency  modes  can  be  captured  accurately.  In  facl,  as  examined  by  Penny  et  al.,  bofh  fhe 
Henshell  and  Ong  mefhod  and  fhe  EIDV  mefhod  produce  accepfable  selecfions  in  mosf  cases,  in  a  sub- 
opfimal  manner.  The  DOE  selecfion  algorifhm  based  on  fhe  EIDV  mefhod  is  shown  as  Algorifhm  2.  Eirsf, 
fhe  eigenvalue  problem  Eq.  (16)  is  solved  for  fhe  firsf  k  modes,  and  fhe  associafed  modal  mafrix  is  designafed 
as  =  (<^i,  02,  •  •  • ,  0fc),  or  =  M^fcA^  where  A^  =  diag^^i  The  Eisher  information 

mafrix  A  is  fhen  calculafed  as  A  =  and  an  idempofenf  mafrix  E  is  compufed  as  E  = 

fhe  diagonal  of  which  is  called  fhe  independence  disfribufion  veefor  (see  Ref.^^  for  defailed  formulations.) 
The  leasf  confribufing  DOE  fo  fhe  independence  of  fhe  modes  among  fhe  ones  in  72  is  idenfified  as  fhe  one 
wifh  fhe  smallesf  diagonal  elemenf  in  E.  The  associafed  DOE  are  also  idenfified  and  stored  in  C,  and  bofh 
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Algorithm  2  DOF  selection  based  on  EIDV  method 
1:  Calculate 

2:  for  i  =  1  to  i  =  n7^  —  do 
3: 

4:  E  ^ 

5:  Find  qi  such  that  =  min  e,-,- 

jell 

6:  £  <—  {gi, . . . ,  gnj,}  where  q’s  are  the  DOF  associated  with  the  fcth  contact  pair  (k  €  Cep)  and  is 

the  number  of  DOF  in  /cth  contact  pair 
7:  n\c 

8: 

9:  Delete  rows  of  corresponding  to  the  DOF  in  C 

10:  end  for 


M  and  TZ  are  updated  as  in  the  Henshell  and  Ong  method.  Finally  the  rows  of  corresponding  to  the  DOF 
in  C  are  deleted.  The  iteration  continues  until  the  size  of  M  reaches  the  desired  number  nj\f. 

Although  the  EIDV  method  shares  similar  objective  for  choosing  DOE  with  the  Henshell  and  Ong 
method,  the  objective  of  the  EIDV  method  is  not  exactly  the  maximization  problem  of  Eq.  (18).  Therefore 
it  is  expected  that  the  Henshell  and  Ong  method  returns  better  solutions  to  the  given  maximization  problem 
than  the  EIDV  method,  as  it  is  shown  in  the  next  section. 

IV.  Case  studies 

In  section  III,  the  method  to  select  the  nonlinear  DOE  has  been  introduced.  In  this  section,  the  validity 
and  applicability  of  the  method  are  discussed  by  applying  the  algorithm  to  two  example  problems.  With 
the  first  case  study,  the  validity  of  the  proposed  method  is  discussed  in  terms  of  the  bilinear  frequencies  and 
forced  response.  Eurthermore  a  metric  to  assess  their  accuracy  is  introduced  and  examined.  The  second  case 
study  is  provided  to  demostrate  the  applicability  of  the  proposed  method  to  the  system  with  a  more  realistic 
EE  model  featuring  a  large  number  of  DOE  on  the  faces  involving  intermittent  contact. 

IV.A,  Simple  cracked  plate  model 

IV.A.l.  Problem  description 

A  cantilevered  cracked  plate  model  was  constructed  with  Young’s  modulus  E  =  2.0  x  10^^  Pa,  Poisson’s 
ratio  u  =  0.3,  and  density  p  =  7800  kg/m^,  and  its  geometry  is  shown  as  Pig.  3a  where  w  =  6.0  x  10“^m, 
I  =  6.0  X  10“^m,  h  =  1.5  X  10“^m,  Ic/l  =  0.625,  and  hc/h  =  0.475.  The  model  was  discretized 
with  5,120  linear  solid  elements  and  resulted  in  mass  and  stiffness  matrices  with  18,630  DOE  On  the  crack 
surfaces  as  shown  in  Pig.  3b,  there  are  180  nodes,  or  90  contact  pairs  on  the  surfaces  hence  the  number  of 
the  associated  DOE  is  540.  The  CMS  method  shown  in  the  section  UFA  was  then  applied  to  the  PE  model, 
and  it  resulted  in  the  681  DOE  (3.6%  of  the  original  size)  system  consisting  of  621  physical  DOE  and  60 
modal  coordinates  corresponding  to  the  free-interface  normal  modes.  With  this  reduced  order  model,  both 
algorithms  in  Algorithms  1  and  2  were  applied  for  nj^  =4,  8,  16,  32,  64,  and  128.  Por  the  EIDV  algorithm, 
the  first  four  modes  were  considered  to  construct  the  modal  matrix,  or  $4. 

In  order  to  compare  these  results  with  an  “intuitive”  selection  method,  a  selection  criteria  was  also 
employed,  where  the  nonlinear  DOE  were  chosen  based  on  the  amount  of  penetration  between  the  nodes  in 
a  contact  pair  for  the  modes  of  interest,  which  in  this  case  is  the  fourth  mode.  Namely,  it  was  hoped  that 
penalizing  the  inter-penetration  of  the  most  penetrating  contact  pairs  would  produce  the  “stiffest”  system 
response.  The  selected  node  pattern  with  such  criterion  is  shown  as  Pig.  4,  and  the  results  of  the  EIDV 


11  of  21 


American  Institute  of  Aeronautics  and  Astronautics 


Figure  3.  Cantilevered  cracked  plate  model:  (a)  FE  model,  (b)  Magnified  crack  surface 


(a)  nAA  =  4  (2  pairs) 


(b)  nAT  =  8  (4  pairs) 


(c)  nj\f  =  16  (8  pairs) 


(d)  nj\f  =  32  (16  pairs)  (e)  nj^  =  64  (32  pairs)  (f)  nj\f  =  128  (64  pairs) 

Figure  4.  Selected  nodes  by  an  intuitive  approach  (left  edge  open) 


method  and  the  Henshell  and  Ong  method  are  shown  as  Figs.  5  and  6.  As  can  be  seen  in  Fig.  4,  if  the  nodes 
are  chosen  based  on  the  amount  of  penetration,  the  selection  starts  from  the  nodes  near  the  crack  edge  (open 
side)  for  n_sf  =  4,  and  it  then  proceeds  toward  the  tip  of  the  crack  (closed  side)  as  nj\f  increases.  It  makes 
sense  because  the  motion  of  the  crack  surfaces  is  significant  near  the  open  edge  than  that  near  the  closed 
edge.  On  the  other  hand  with  the  EIDV  method,  the  method  also  starts  to  select  the  nodes  near  the  crack 
edge,  but  it  tends  to  choose  more  nodes  on  the  crack  rims  than  the  nodes  near  the  crack  edge  as  shown  in 
Fig.  5.  Finally  with  the  Henshell  and  Ong  method,  it  also  select  the  nodes  near  the  crack  edge  first,  for 
nj\f  =  4,  but  it  then  tends  to  select  the  nodes  over  the  crack  surface  in  a  more  distributed  manner  as  can  be 
seen  in  Fig.  6. 

IV.A.2.  Forced  Response  Analysis 

Next,  in  order  to  see  the  influence  of  fhe  application  of  fhe  nonlinear  B.C.  onfo  fhe  selecfed  nodes  on  an 
NNM  frequency,  forced  response  analysis  was  carried  ouf  by  applying  an  exfernal  harmonic  loading  fo  fhe 
cracked  plate.  As  one  mighf  nofice,  when  fhe  forced  response  of  fhis  sfrucfure  wifh  a  crack  is  considered,  fhe 
repefifive  opening  and  closing  of  fhe  crack  faces  musf  be  freafed  appropriafely  wifh  confacf  algorifhms.  As  a 
resulf,  the  vibration  is  nonlinear  and  the  steady-state  response  of  the  displacement  may  not  be  expressed  as  a 
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(a)  tiaa  =  4  (2  pairs) 


(d)  TiAf  =  32  (16  pairs) 


(b)  nAT  =  8  (4  pairs) 


(e)  n/yf  =  64  (32  pairs) 


(c)  n^f  =  16  (8  pairs) 


(f)  n^f  —  128  (64  pairs) 


Figure  5.  Selected  nodes  by  EIDV  method  (left  edge  open) 


(a)  nj</  =  4  (2  pairs)  (b)  nj\r  =  8  (4  pairs)  (c)  =  16  (8  pairs) 


(d)  nj\f  =  32  (16  pairs)  (e)  nj^  =  64  (32  pairs)  (f)  nj^  =  128  (64  pairs) 

Figure  6.  Selected  nodes  by  the  modified  Henshell  and  Ong  method  (left  edge  open) 


harmonic  function  even  if  the  external  force  is  a  harmonic  function.  Therefore  in  this  study,  the  steady-state 
response  was  obtained  by  assuming  that  the  displacement  can  be  expressed  as  a  truncated  Fourier  series,  and 
the  nonlinear  boundary  condition  can  be  enforced  by  the  penalty  method.^  The  method  is  called  the  hybrid 
frequency-time  domain  method,^®’ which  is  based  on  the  concept  of  harmonic  balance  method."*  Detailed 
formulation  of  the  method  is  omitted  in  this  paper. 

It  is  noted  that  the  system  matrices  were  further  reduced  by  the  application  of  Eq.  (4)  to  the  reduced- 
order  model  before  the  forced  response  calculation,  by  keeping  the  selected  node  pairs  as  the  active  DOF 
and  condensing  out  the  other  DOF  including  physical  and  modal  coordinates.  For  example,  with  =  64 
(32  pairs),  the  system  size  was  reduced  down  to  155  DOF,  which  is  0.83%  of  the  original  system  size. 

A  harmonic  forcing  of  magnitude  3  N  was  then  applied  at  the  tip  of  the  plate,  in  order  to  excite  the 
first  vibration  mode,  which  corresponds  to  the  first  out-of-plane  bending  mode.  The  forced  response  was 
calculated  for  both  linear  case,  i.e.,  with  the  open  B.C.,  and  nonlinear  case  with  the  nonlinear  boundary 
conditions  imposed  on  M  with  the  selections  by  the  Henshell  and  Ong  method.  The  results  are  shown  in 
Figs.  7a  and  7b.  As  can  be  seen  in  Fig  7a,  the  selection  pattern  does  not  alter  the  linear  forced  response. 
This  is  because  the  selection  of  the  active  DOF  does  not  alter  the  eigenvalues  of  the  reduced  order  model. 


Table  1.  Average  CPU  time  in  second  per  frequency  to  obtain  forced  response  for  256  steps  witbin  200  to  215  Hz 


nM 

Intuitive 

EIDV 

Henshell  &  Ong 

2 

4.98x10“** 

4.63x10“** 

4.73x10“** 

4 

2.21x10-2 

1.57x10-2 

1.27x10-2 

8 

7.27x10-2 

1.13x10“* 

8.13x10-2 

16 

1.96x10“* 

4.19x10“* 

3.67x10“* 

32 

5.79 

5.38 

3.20 

64 

4.24x10 

3.87x10 

4.04x10 

90  (full) 

1.02  X  102 
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Figure  7.  Results  of  forced  response  analysis  of  the  cracked  plate 


and  it  was  assumed  that  the  system  was  completely  linear  when  the  forced  response  was  calculated.  On  the 
other  hand,  the  number  of  contact  pairs  greatly  affects  the  results  of  nonlinear  forced  response  as  shown  in 
Fig.  7b.  One  may  observe  that  the  response  with  64  contact  pairs  is  almost  identical  to  that  with  the  full 
set  of  90  contact  pairs,  which  implies  that  for  accurately  calculating  the  nonlinear  resonant  frequencies,  it 
may  not  be  necessary  to  enforce  nonlinear  boundary  conditions  for  all  the  contact  pairs  on  the  crack  faces. 
The  same  forced  response  calculations  were  carried  out  with  the  node  patterns  selected  by  the  EIDV  method 
and  the  method  based  on  the  amount  of  penetration,  and  they  are  respectively  shown  in  Figs.  7c  and  7d.  As 
can  be  seen  in  Fig.  7c,  the  results  with  the  patterns  chosen  by  the  EIDV  method  are  comparable  with  the 
ones  produced  by  the  Henshell  and  Ong  method.  On  the  other  hand  as  can  be  seen  in  Eig.  7d,  the  forced 
response  with  the  node  patterns  chosen  by  the  “intuitive”  approach  produced  worse  results  than  the  other 
two  methods,  i.e.,  for  a  given  number  of  nj^,  the  predicted  resonant  frequency  by  the  approach  is  lower  than 
that  calculated  by  the  other  methods.  This  is  the  most  visible  in  the  results  for  nj\j-  =  64,  for  which  both  the 
Henshell  and  Ong  method  and  the  EIDV  methods  produced  results  that  are  almost  identical  to  the  results  for 
nj\f  =  90. 

Moreover,  in  order  to  see  the  effects  of  nj^  on  the  computational  speed,  the  average  CPU  times  per  fre¬ 
quency  required  to  obtain  the  steady-state  forced  response  per  frequency,  for  the  frequency  range  considered 
in  Eig. 7  are  shown  in  Table  1.  The  forced  response  was  calculated  at  evenly-spaced  256  points  from  200  to 
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Figure  8.  Errors  in  the  first  four  bilinear  frequencies: (a)  NNM  1,  (b)  NNM  2,  (c)  NNM  3,  (d)  NNM  4 


215  Hz  for  nj^  =  2,  4,  8,  16,  32,  and  64  by  the  three  methods  considered,  and  for  the  case  with  =  90. 
The  CPU  time  was  measured  on  a  Macintosh  computer  with  Intel  Core  2  Duo  2.4GHz  and  4.0GB  of  RAM. 
Although  there  is  no  significant  difference  between  the  selection  algorithms,  we  can  see  that  the  CPU  time 
can  be  greatly  reduced  when  the  calculations  were  done  with  the  sampled  nodes,  without  sacrificing  the 
accuracy  if  the  proper  set  of  nodes  were  chosen  by  the  proposed  method.  It  is  noted  that  the  average  CPU 
time  is  approximately  a  cubic  function  of  nj^. 

IV. A.  3.  Bilinear  Frequency  Approximation 

Finally,  the  influence  of  the  selected  node  pattern  on  the  bilinear  frequencies  is  discussed.  The  first  four 
bilinear  frequencies  were  calculated  for  the  model  with  the  selected  node  patterns  with  the  three  methods, 
and  the  results  are  shown  in  Fig.  8  along  with  their  linear  vibration  mode  shapes.  The  first  four  modes 
correspond  to  the  first  out-of-plane  bending,  the  first  torsion,  the  second  out-of-plane  bending,  and  the  first 
in-plane  bending  modes  respectively.  The  plots  in  Fig.  8  show  the  percentage  errors  in  the  bilinear  frequency 
versus  the  number  of  contact  pairs,  where  the  error  is  defined  as  the  ratio  of  the  difference  between  the 
bilinear  frequency  with  the  sampled  contact  pairs  and  that  with  the  full  set  of  contact  pairs,  to  that  with  the 
full  set  of  contact  pairs.  As  can  be  seen,  the  Henshell  and  Ong  method  consistently  provides  the  best  results 
among  all  the  methods  for  the  first  four  modes.  Moreover,  it  shows  the  best  convergence  rate  in  terms  of  the 
number  of  contact  pairs. 

IV.A.4.  A  posteriori  accuracy  assessment 

As  seen  above,  even  though  the  intuitive  approach  chooses  the  contact  pairs  that  show  the  most  penetration, 
application  of  the  nonlinear  boundary  conditions  to  these  nodes  does  not  result  in  the  “stiffest”  vibration 
response.  To  be  specific,  the  Henshell  and  Ong  method  and  the  EIDV  method  produced  the  node  patterns 
that  yield  the  closer  results  to  the  reference  results  in  terms  of  forced  response  and  bilinear  frequencies,  than 
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Figure  9.  Virtual  impulse  for  a  period  of  vibration  for  NNM  1 


the  patterns  chosen  by  the  intuitive  approach.  In  particular,  the  Henshell  and  Ong  method  iteratively  aims 
to  solve  the  maximization  problem  Eq.  ( 1 8)  in  a  sub-optimal  manner.  Therefore  the  bilinear  frequencies  as 
well  as  the  resonant  frequencies  were  well  approximated  with  the  nodes  chosen  by  the  Henshell  and  Ong 
method.  In  order  to  better  understand  the  governing  factor  for  the  accuracy  of  the  results,  a  more  physical 
interpretation  of  the  results  is  provided  here.  Namely,  the  key  effect  for  achieving  the  good  approximation 
of  the  NNM  frequency  is  to  ensure,  as  much  as  possible,  the  non-penetrability  condition  on  the  contact  pairs 
where  the  nonlinear  B.C.s’  are  not  applied.  The  penetration  should  be  evaluated  during  a  vibration  cycle, 
thus  both  the  depth  and  the  duration  of  the  penetration  should  be  taken  into  account.  These  quantities  vary  in 
space,  and  depend  on  the  frequency  of  vibration.  Hence  as  a  metric  to  characterize  not  only  the  amount  but 
also  the  duration  of  penetration  over  the  entire  crack  surfaces  for  a  given  vibration  frequency,  the  following 
quantity  is  introduced: 

F  =  f  (  [  A:ettp(r,f)dr  )  dt  (23) 

Jo  \JrA(rB)  J 

where  F  is  a  quantity  with  the  dimension  of  impulse  named  “virtual  impulse”,  kg  is  an  equivalent  spring 
constant  per  unit  length  determined  by  the  ratio  between  the  Young’s  modulus  multiplied  by  the  characteris¬ 
tic  area  and  the  characteristic  length,  Up  is  the  amount  of  penetration  along  the  surface  normals,  and  T  is  the 
period  of  vibration.  The  quantity  F  is  calculated  based  on  the  calculated  time  trajectory  of  displacements  of 
the  nodes  on  the  crack  surfaces,  and  can  be  thought  of  as  an  impulse  that  does  not  contribute  to  the  system 
response,  as  this  impulse  is  not  applied  to  the  system  when  the  response  is  calculated.  In  other  words,  the 
smaller  the  value  of  F  is,  the  stricter  the  boundary  conditions  are  imposed  on  the  nodes  over  the  entire  crack 
surfaces. 

First,  the  forced  response  analysis  was  carried  out,  and  the  corresponding  time  history  of  Up  over  the 
entire  crack  surface  was  recovered  from  the  vibration  response.  The  integrals  in  the  Eq.  (23)  were  then 
evaluated  by  a  simple  quadrature  rule  both  in  space  and  time.  The  metric  was  calculated  for  the  first  and 
the  fourth  modes  for  8  and  32  pairs  chosen  by  the  methods,  and  the  results  are  shown  in  Figs.  9  and  10. 
As  can  be  seen  in  Figs.  9  and  10,  the  virtual  impulse  varies  over  the  frequency  range.  In  particular,  when 
the  frequency  of  excitation  is  close  to  the  resonant  frequency,  or  the  NNM  frequency,  then  the  amount  of 
penetration  increases  as  well.  However  for  all  cases,  the  Henshell  and  Ong  method  consistently  results  in  the 
smallest  impulse  over  the  frequency  range  among  the  three  methods  considered.  It  means  that  the  nonlinear 
B.C.  on  the  crack  faces  is  the  most  strictly  enforced  by  the  node  patterns  chosen  by  the  Henshell  and  Ong 
method. 
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IV.B.  Large  scale  cracked  blade  model 

In  this  subsection,  the  application  to  the  proposed  method  to  a  finite  element  model  of  a  cracked  blade  model 
with  a  large  number  of  DOF  is  considered.  A  blade  model  whose  thickness  is  2.5  x  10“^m  and  both  chord 
and  span  lengths  are  approximately  5.0  x  10~^m  was  constructed,  with  the  Young’s  modulus  E  =  205GPa, 
density  p  =  7832kg/m^,  and  Poisson’s  ratio  u  =  0.3.  And  it  was  discretized  with  linear  and  quadratic 
tetrahedral  elements,  resulting  in  a  system  with  392,361DOF  as  shown  in  Fig.  11a.  The  crack  path  as  well 
as  its  surrounding  finite  element  mesh  was  generated  by  a  fracture  analysis  code  FRANC3D.‘^‘^  The  crack 
surfaces  consist  of  quadratic  elements  with  487  contact  pairs  at  the  edges  of  vertices  of  triangles,  resulting 
in  2,922  DOF  for  both  surfaces,  as  shown  in  Fig.  1  lb.  The  size  of  the  FE  model  was  then  reduced  down  to 
2,949  DOF  (0.75%  of  the  original  size)  by  the  primary  CMS  method,  which  consists  of  2,934  physical  DOF 
(2,922  DOF  on  crack  faces  and  additional  12  DOF  for  forcing)  and  15  modal  coordinates  corresponding  to 
the  free-interface  normal  modes.  The  proposed  method  was  then  applied  to  the  reduced  order  model,  and  the 
results  are  shown  in  Figs.  12  and  13.  For  the  EIDV  method,  the  first  15  modes  were  used  for  the  calculations. 
As  can  be  seen  in  Eig.  12,  the  EIDV  method  again  selected  the  nodes  along  the  rim  of  the  crack  faces,  which 
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Figure  12.  Selected  nodes  by  EIDV  method  for  ng  =  974  (487pairs) 


(c)  riAT  =  256  (128  pairs)  (d)  nj\f  =  512  (256  pairs) 

Figure  13.  Selected  nodes  by  the  modified  Henshell  and  Ong  method  for  ng  =  974  (487pairs) 


is  similar  to  the  previous  case  study.  On  the  other  hand,  the  Henshell  and  Ong  method  tends  to  choose 
the  nodes  slightly  off  of  the  crack  rim  in  a  more  distributed  manner  over  the  crack  face,  as  in  the  previous 
example.  With  the  selected  node  patterns,  the  first  four  bilinear  frequencies  were  calculated.  The  first  four 
modes  correspond  to  the  first  out-of-plane  bending,  the  first  torsion,  the  second  out-of-plane  bending,  and 
the  first  chord-wise  bending  mode  respectively.  The  errors  were  then  calculated  as  in  the  previous  example 
and  shown  in  Fig.  14,  along  with  their  linear  vibration  mode  shapes.  As  can  be  seen,  the  Henshell  and  Ong 
method  consistently  shows  the  better  results  than  the  results  obtained  by  the  EIDV  method. 

V.  Conclusion 

In  this  paper,  a  novel  reduced  order  modeling  framework  for  the  nonlinear  vibration  analysis  of  elastic 
structures  with  intermittent  contact  was  proposed.  In  section  III,  the  modeling  framework  was  developed 
based  on  a  method  of  component  mode  synthesis.  The  master  DOF  selection  scheme  for  Guyan  reduction 
was  formulated  by  considering  the  close  relationship  between  the  optimal  master  DOF  selection  and  the  op¬ 
timal  constraint  locations  for  maximizing  the  fundamental  natural  frequency.  The  method  is  a  combination 
of  the  sequential  elimination  method  proposed  by  Henshell  and  Ong,  and  appropriate  coordinate  transfor¬ 
mations  to  the  reduced  order  model.  Another  method  for  choosing  the  nodes  was  also  introduced  for  the 
sake  of  comparative  study,  which  is  based  on  a  method  to  optimally  choose  the  measurement  locations  such 
that  the  measured  mode  becomes  as  linearly  independent  as  possible.  The  method  was  then  applied  to  a 
representative  finite  element  model  in  section  IV.  In  IV.A,  the  methods  were  applied  to  a  cracked  plate 
model.  Using  the  selected  node  patterns,  forced  response  analysis  was  carried  out  to  see  the  effects  of  the 
selection  patterns  on  the  frequency  response.  Furthermore  the  resonant  frequencies  were  calculated  by  the 
application  of  bilinear  frequency  approximation.  It  was  confirmed  that  the  selected  DOF  resulted  in  accu- 
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(a)  (b) 


(c)  (d) 

Figure  14.  Error  in  the  bilinear  frequency  for  (a)  NNM  1,  (b)  NNM  2  (c)  NNM  3  (d)  NNM  4 


rate  prediction  of  nonlinear  resonant  frequencies  in  comparison  to  the  benchmark  case  of  using  all  DOF  on 
the  crack  surfaces.  Furthermore,  it  was  demonstrated  that  the  method  also  achieves  the  reduction  in  CPU 
time  for  the  nonlinear  forced  response  calculations,  without  sacrificing  the  accuracy  in  the  predicted  forced 
response.  Moreover,  a  posteriori  accuracy  assessment  procedure  was  introduced  by  examining  the  amount 
of  penetration  on  the  crack  surfaces  during  a  vibration  cycle.  In  IV.B,  the  method  was  also  applied  to  an 
academic  cracked  blade  model  with  a  much  larger  number  of  DOF  on  crack  faces.  The  node  selection  pat¬ 
terns  as  well  as  the  errors  in  the  bilinear  frequencies  conform  to  the  results  in  the  example  in  IV. A.  For  the 
methods  examined,  the  node  patterns  selected  using  the  proposed  new  method  consistently  showed  the  best 
results. 
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